      SUBROUTINE DEREQ
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 MPRIME,NSUBX,JARRAY,ISUBD
      INTEGER*4 ACNTRL
C
C     SUBROUTINE DEREQ1 IS THE CENTER OF THE FSD PROGRAM. IT SETS
C     UP THE EQUATIONS OF MOTION AND DIRECTS THE CALCULATION OF
C     THE DERIVATIVES OF THE STATE VARIABLES, CALL ROUTINES TO
C     COMPUTE THE CONTROL AND EXTERNAL TORQUES
C
C     IT IS CALL BY  ADMIMP  THE ADAMS-MOULTON FOURTH ORDER
C     PREDICTOR AND CORRECTOR INTEGRATOR
C     DEREQ1 IS ORIGINALLY WRITTEN BY E.A. LAWLOR OF AVCO SYSTEM
C     DIVISION.  MODIFIED BY K. YONG OF COMPUTER SCIENCES CORP.
C     FOR VARIOUS ADDITIONAL CAPABILITIES
C
C
      COMMON/CONDMP/ IDPHLD
C
      COMMON/CONSTS/ PI,TWOPI,RADIAN
C
      COMMON/CRATIO/ RATIO
C
      COMMON/DEBUG2/ IOUT,JOUT,KLUGE
C
      COMMON/DMMNT1/ ZKBM(6),FMAK(10),FMBK(10),ARTETA(3),CMTORK(3),
     .               ITORK,IBENDM,ITENSE,ITNS1
C
      COMMON/GRNTST/ ALFAEG,DELTAG,PHASEG,ALTUDE,OMGY(3),
     *               GACC(3),GLOCAT(3),IGRUND,IALTUD,IGASBR
C
      COMMON/HVCOMP/ YBAR(3),YBARD(3),CIN(3,3),CID,CIND(3)
C
      COMMON/IMAIN1/ IDATE,LSAVE,INOPT,IPLOT,NUMEQS,IPLTPE,IORB,ITAPE
C
      COMMON/INEWR / NKT(10),ICP,ICPS
C
      COMMON/INTTRP/ ITPROT,NUMTIP(10),IRDBUG
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      COMMON/ITCNTL/ IPULSE,ISPLSE,KPULSE,ITSW,IOTSW
C
      COMMON/ITW   / ITWIST,ITWST1
C
C
      COMMON/MOMENT/ ACNTRL,IVISCS,IATTDE,IMGMTS,IWHEEL,NPULSE
C
      COMMON/NODER / NDER,NOPT
C
      COMMON/OUTFOR/ SUMMTS(3)
C
      COMMON/RATTDE/ DTMXA,PXI,PXO,EIT,REFTIM,TMX1,TMX2,CMX,ISW
C
      COMMON/RMAIN1/ DELTAT,FACTOR,FREQ,TSTOP,DELMIT,
     .               UPBND(150),DNBND(150)
C
      COMMON/RPOOL1/ RHOK(10),TIM1,SA(3,3),FM1(3,3),ZLK(10),OMEG(3),
     .               ZLKP(10),ZLKDP(10),CMAT(3,3),GBAR(3,3),YBCM(3),
     .               ZBZK(3,10),FCM(3,3),DTO,PHID,PHI
C
      COMMON/RPOOL2/ PO,SD(3),DAN(3,10),DBN(3,10),CFMT(3,3),DIY1(3),
     .               SD1(3),DT1,P1,AERO,DTO1,YIZK(3),PO1
C
      COMMON/RPOOL9/ RT1(7),RT2(10,9),ALP(7,7),GAM(10,9,7),DEL(10,9,9)
C
      COMMON/RVISCS/ NSUBX(3),YARRAY(3),RADTBE,VISCTY,RADRNG,DENSTY,
     .               JARRAY(3),SSUBY,OMEGL,ISUBD
C
      COMMON/SATPOS/ TLAST,TFRST,XLAST(3),XFRST(3),XDLAST(3),XDFRST(3),
     .               SDLAST(3),SDFRST(3),ADLAST(3),ADFRST(3),XDT
C
      COMMON/TRQOUT/ OUTTRQ(150)
C
      COMMON/VARBLS/ DEPEND(150),DERIV(150)
C
      COMMON/VECTRS/ XSAT(3),XSATDT(3),AD(3)
C
      COMMON/ZSPINR/ DTZMA,PZDT,CMZO,ISPIN3,JSPIN
C
C
      DIMENSION DDEL(9,9),DELINV(10,9,9),ETA(7),ZML(7,7),RT(9)
C
      DIMENSION C(10),DUM(7,9),MT(10),VPRIME(3,2)
C
      DIMENSION OMEGDT(3),RWHEEL(3),XMB(3),MPRIME(3),YDOT(3)
C
      DIMENSION GAME(9,7)
C
      DIMENSION RVEC(28),AMAT(28),AUX(7)
C
      REAL*4 EPS/1.0E-6/
C
C        SPECIAL DEBUGGING OUTPUT
C
      NDER=NDER+1
      IF(NOPT.EQ.1) WRITE(6,302) NDER
  302 FORMAT(15X,' NDEREQ1 = ',I5)
      IF(JOUT.NE.2) GO TO 4
      IOUT=1
      TOUT=DMOD(TIM1,8.64D4)
      IF(TOUT.GT.FREQ) TOUT=DMOD(TOUT,FREQ)
      IF(TOUT.EQ.FREQ.OR.(FREQ-TOUT).LE.DELTAT) IOUT=2
    4 IF(JOUT.EQ.3) IOUT=2
C
C        SETUP FOR CONTROL FORCES,WHEELS,ATTITUDE CONTROL,VISCOUS
C        NUTATION DAMPING,SPIN UP TORQUEING
C
      N=NUMEQS - 4
      CMXO=0.D0
      DO 5 I=1,3
      RWHEEL(I)=0.D0
      XMB(I)=0.D0
      MPRIME(I)=0.D0
      VPRIME(I,1)=0.0D0
      VPRIME(I,2)=0.0D0
      IF(IVISCS.EQ.0) GO TO 5
      LL=N+I
      YARRAY(I)=DEPEND(LL)
    5 CONTINUE
      IF(IVISCS.EQ.1) SSUBY=DEPEND(NUMEQS)
      CMZO1=0.D0
      IF(ISPIN3.NE.0.AND.JSPIN.NE.0) CMZO1=CMZO
C
      DO 15 I=1,7
      DO 14 J=1,7
      ZML(I,J)=0.0D0
   14 CONTINUE
      ETA(I)=0.0D0
   15 CONTINUE
C
      IF(IOUT.EQ.1) GO TO 16
      WRITE(6,7001) ZLK
      WRITE(6,7002) ZLKP
 7001 FORMAT('0',1X,'ZLK(K)  1 ',1P10E12.4)
 7002 FORMAT('0',1X,'ZLKP(K) 1 ',1P10E12.4)
C
   16 CONTINUE
C
C        OBTAIN ORBITAL DATA
C
      CALL SETVAL(2)
C
C        ITIM IS SWITCH FOR CALCULATING - NEEDED AT TIME ZERO ONLY
C
      IF(ITIM.EQ.2) GO TO 20
      CALL COMBNZ
      GO TO 40
C
C        IF LENGTH IS CHANGING AFTER  TIME ZERO COMBNZ IS CALLED AT
C        EACH TIME
C
   20 DO 30 K=1,IK
      IF(ZLKP(K).EQ.0.D0) GO TO 30
      CALL COMBNZ
      GO TO 40
   30 CONTINUE
C
   40 IF(IOUT.EQ.1) GO TO 17
      WRITE(6,7003) ZLK
      WRITE(6,7004) ZLKP
 7003 FORMAT('0',1X,'ZLK(K)  2 ',1P10E12.4)
 7004 FORMAT('0',1X,'ZLKP(K) 2 ',1P10E12.4)
C
   17 CONTINUE
C
      CALL FINDGB
C
      CALL THERME
C
      CALL FNDALP
C
      IF(ITWIST.NE.0.OR.ITPROT.NE.0) CALL DEREQT(DEPEND)
C
      DO 35 I=1,3
      I3=I+3
      DO 35 J=1,3
      J3=J+3
      CIN(I,J)=ALP(I3,J3)
   35 CONTINUE
C
      IF(IOUT.EQ.1) GO TO 41
C
C        DEBUGGING OUTPUT FOR IOUT=2
C                       TEMPORARY WRITE'S
      WRITE(6,7005) ZLK
      WRITE(6,7006) ZLKP
 7005 FORMAT('0',1X,'ZLK(K)  3 ',1P10E12.4)
 7006 FORMAT('0',1X,'ZLKP(K) 3 ',1P10E12.4)
C
      T=DMOD(TIM1,8.64D4)
      WRITE(6,8000)
 8000 FORMAT(15X,'DEBUGGING OUTPUT FROM DEREQ1')
      WRITE(6,8001)
 8001 FORMAT(15X,' SA MATRIX')
      WRITE(6,9000) ((SA(I,J),J=1,3),I=1,3)
 9000 FORMAT(' ',2X,1P9E14.6//)
      WRITE(6,9001) (RT1(I),I=1,7)
 9001 FORMAT(' ',2X,'RT1',7E16.8/)
      WRITE(6,8002)
 8002 FORMAT(15X,' RT2 ARRAY ')
      WRITE(6,9002) ((RT2(I,J),J=1,9),I=1,10)
 9002 FORMAT(9E14.6)
      WRITE(6,8003)
 8003 FORMAT(15X,'DEPEND ARRAY FROM DEREQ1')
      WRITE(6,9002) (DEPEND(I),I=1,150)
      WRITE(6,8004)
 8004 FORMAT(15X,'DERIV ARRAY FROM DEREQ1')
      WRITE(6,9002) (DERIV(I),I=1,150)
      WRITE(6,9003) T,LSAVE,DELTAT
 9003 FORMAT('0',2X,'DEREQ - T=',E20.8,2X,'LLL=',I2,2X,'DELTAT=',E20.8)
      WRITE(6,8006)
 8006 FORMAT(15X,' FM1 ARRAY ')
      WRITE(6,9000) ((FM1(I,J),J=1,3),I=1,3)
C
   41 CONTINUE
      DO 50 I=1,7
      ETA(I)=RT1(I)
      DO 50 J=1,7
   50 ZML(I,J)=ALP(I,J)
C
      MAXD=9
      NALP=6
      IF(IDAMP.EQ.1) NALP=7
C
      DO 42 I=1,28
      RVEC(I)=0.0D0
   42 CONTINUE
C
      NRHS=1
      I1=NALP
      IF(IDPHLD.EQ.0) GO TO 43
      NRHS=NRHS+1
      I1=I1+NALP
      RVEC(I1)=1.0D0
   43 CONTINUE
      IF(IGRUND.EQ.0.OR.IGASBR.EQ.0) GO TO 44
      NRHS=NRHS+2
      I1AX=I1+4
      I2AX=I1AX+NALP+1
      RVEC(I1AX)=1.0D0
      RVEC(I2AX)=1.0D0
   44 CONTINUE
C
      CALL ACTFLT(1,TIM1,ETA)
C
      DO 130 K=1,IK
C
      CON=RHOK(K)*ZLK(K)
C
      NKB=NK(K)
      NTW=NKT(K)
      NSZ=2*NKB+NTW
C
      IF(NSZ.EQ.0) GO TO 130
C
      DO 60 I=1,NSZ
      DO 58 J=1,NSZ
      DDEL(I,J)=DEL(K,I,J)
   58 CONTINUE
      DO 59 J=1,7
      GAME(I,J)=GAM(K,I,J)
   59 CONTINUE
   60 CONTINUE
C
      IF(IOUT.EQ.1) GO TO 62
C
      WRITE(6,6010) ((DDEL(I,J),J=1,9),I=1,NSZ)
 6010 FORMAT(' DELTA         ',1P9E13.6)
      WRITE(6,6011) ((GAME(I,J),J=1,7),I=1,NSZ)
 6011 FORMAT(' GAMMA         ',1P7E13.6)
C
   62 CONTINUE
C
C     INVERSE OF DELTA MATRIX
C
      CALL INVERT(DDEL,NSZ,MAXD,IND,C,MT)
C
C
C
      IF(IOUT.EQ.1) GO TO 64
      WRITE(6,6012) ((DDEL(I,J),J=1,9),I=1,NSZ)
 6012 FORMAT(' DELTA INVERSE ',1P9E13.6)
C
   64 CONTINUE
C
C     STORE DELTA INVERSE
C
      DO 70 I=1,NSZ
      DO 70 J=1,NSZ
      DELINV(K,I,J)=DDEL(I,J)
   70 CONTINUE
C
      DO 80 I=1,NALP
      DO 80 J=1,NSZ
      DUM(I,J)=0.D0
      DO 80 L=1,NSZ
   80 DUM(I,J)=DUM(I,J) + GAM(K,L,I)*DDEL(L,J)
C
C
      DO 100 I=1,NALP
      DO 100 J=1,NALP
      SUM1=0.D0
      DO 90 L=1,NSZ
   90 SUM1=SUM1 + DUM(I,L)*GAM(K,L,J)
  100 ZML(I,J)=ZML(I,J) - SUM1*CON
C
C
      DO 120 I=1,NALP
      SUM2=0.D0
      DO 110 L=1,NSZ
  110 SUM2=SUM2 + DUM(I,L)*RT2(K,L)
  120 ETA(I)=ETA(I) - SUM2*CON
  130 CONTINUE
C
C
      IF(KPULSE.NE.IPULSE) CALL PULSER(ETA)
C
C     CALL GRDFCE TO OBTAIN THE GRAVITATIONAL FORCE ON THE GROUND
C
      IF(IGRUND.EQ.1) CALL GRDFCE(ETA)
C
C
      IF(IOUT.EQ.1)GO TO 131
      WRITE(6,7007) ZLK
      WRITE(6,7008) ZLKP
 7007 FORMAT('0',1X,'ZLK(K)  4 ',1P10E12.4)
 7008 FORMAT('0',1X,'ZLKP(K) 4 ',1P10E12.4)
C
      WRITE(6,9500) (ETA(I),I=1,7)
      WRITE(6,9501)
      WRITE(6,9502) ((ZML(I,J),J=1,7),I=1,7)
 9500 FORMAT('0','ETA=',7E17.8)
 9501 FORMAT('0','ZML')
 9502 FORMAT(' ',5X,7E17.8)
      WRITE(6,8005)
 8005 FORMAT(15X,' END OF DEBUGGING OUTPUT FROM DEREQ1')
C
  131 CONTINUE
C
C
C                       TEST FOR ADDING MOMENTUM WHEELS
C
      IF(IWHEEL.EQ.1) CALL WHEELS(1,RWHEEL)
C
      CALL WHREAC(ETA)
      CALL AWREAC(ETA)
C
      OUTTRQ(1)=RWHEEL(1)
      OUTTRQ(2)=RWHEEL(2)
      OUTTRQ(3)=RWHEEL(3)
C
C                       TEST FOR MAGNETIC MOMENTS
C
      IF(IMGMTS.NE.0) CALL MAGNTS(XMB,SA)
C
C                       TEST FOR ATTITUDE CONTROL
C
      IF(NPULSE.EQ.0) GO TO 140
      IF(IATTDE.EQ.1) CALL XACM(CMXO)
C
C
C                       TEST FOR VISCOUS DAMPER
C
  140 CONTINUE
      IF(IVISCS.EQ.1) CALL VISCUS(1,ZML,OMEG,OMEGDT,SYDOT,YDOT,MPRIME)
C
      CALL VISCS2(1,ZML,OMEG,OMEGDT,DEPEND,DERIV,VPRIME)
C
C
C                       ADD ATTITUDE CONTROL MOMENT TO X-AXIS(ETA(4))
C
      ETA(4)=ETA(4) + CMXO
      ETA(6)=ETA(6) + CMZO1
C
C                       ADD OTHER MOMENTS TO ETA
C
      ETA(4)=ETA(4) + VPRIME(1,1) + VPRIME(1,2)
      ETA(5)=ETA(5) + VPRIME(2,1) + VPRIME(2,2)
      ETA(6)=ETA(6) + VPRIME(3,1) + VPRIME(3,2)
C
      ETA(4)=ETA(4) + RWHEEL(1) + XMB(1) + MPRIME(1)
      ETA(5)=ETA(5) + RWHEEL(2) + XMB(2) + MPRIME(2)
      ETA(6)=ETA(6) + RWHEEL(3) + XMB(3) + MPRIME(3)
C
C                       SAVE THE TOTAL MOMENTS FOR OUTPUT
C
      SUMMTS(1)=SUMMTS(1)+XMB(1)+CMXO
      SUMMTS(2)=SUMMTS(2)+XMB(2)
      SUMMTS(3)=SUMMTS(3)+XMB(3)+CMZO1
C
      DO 150 I=1,NALP
      IF(ZML(I,I).EQ.0.D0) GO TO 150
      DO 150 J=1,NALP
      IF(DABS(ZML(I,J)/ZML(I,I)).LT.RATIO) ZML(I,J)=0.D0
  150 CONTINUE
C
      IF(ITORK.EQ.0) GO TO 759
      ETA(4)=ETA(4)+CMTORK(1)
      ETA(5)=ETA(5)+CMTORK(2)
      ETA(6)=ETA(6)+CMTORK(3)
      SUMMTS(1)=SUMMTS(1)+CMTORK(1)
      SUMMTS(2)=SUMMTS(2)+CMTORK(2)
      SUMMTS(3)=SUMMTS(3)+CMTORK(3)
  759 CONTINUE
C
      CALL SECBD2(1,ZML,ETA,NALP)
C
      CALL PCSNSR
C
      CALL DCSNSR
C
      CALL SACSNR
C
      CALL GIMBL2(1,ZML,ETA,NALP)
C
      CALL GMBLD2(1,ZML,ETA)
C
      CALL SAGIM2(1,ZML,ETA,NALP)
C
C
      K=0
      DO 152 I=1,NALP
      DO 151 J=1,I
      K=K+1
      AMAT(K)=ZML(J,I)
  151 CONTINUE
      RVEC(I)=ETA(I)
  152 CONTINUE
C
C
      CALL DGELS(RVEC,AMAT,NALP,NRHS,EPS,IER,AUX)
C
      IF(IER.GE.0) GO TO 153
      WRITE(6,8010) IER
 8010 FORMAT('0 ERROR RETURN FROM DGELS. IER=',I4)
      CALL EXIT
C
  153 CONTINUE
      IF(NRHS.LT.3) GO TO 156
      I1=I1AX
      I2=NALP+I1
      C11=RVEC(I1)
      C21=RVEC(I1+1)
      C12=RVEC(I2)
      C22=RVEC(I2+1)
      R1=RVEC(4)
      R2=RVEC(5)
      DEN=C11*C22-C12*C21
      IF(DEN.EQ.0.0D0) GO TO 155
      GIM1=(C12*R2-C22*R1)/DEN
      GIM2=(C21*R1-C11*R2)/DEN
      DO 154 I=1,NALP
      I1=I1AX-4+I
      I2=NALP+I1
      RVEC(I)=RVEC(I)+RVEC(I1)*GIM1+RVEC(I2)*GIM2
  154 CONTINUE
  155 CONTINUE
      IF(RVEC(4).NE.0.0D0) RVEC(4)=0.0D0
      IF(RVEC(5).NE.0.0D0) RVEC(5)=0.0D0
  156 CONTINUE
C
      IF(IDPHLD.EQ.0) GO TO 702
      IDM=NALP+NALP
      DC=-RVEC(NALP)/RVEC(IDM)
      DO 701 I=1,NALP
      I1=NALP+I
      RVEC(I)=RVEC(I)+DC*RVEC(I1)
  701 CONTINUE
  702 CONTINUE
C
      DO 157 I=1,NALP
      ETA(I)=RVEC(I)
  157 CONTINUE
C
      ARTETA(1)=ETA(1)
      ARTETA(2)=ETA(2)
      ARTETA(3)=ETA(3)
C
C
       DERIV(7)=ETA(4)
       DERIV(8)=ETA(5)
       DERIV(9)=ETA(6)
C
      OMEGDT(1)=ETA(4)
      OMEGDT(2)=ETA(5)
      OMEGDT(3)=ETA(6)
C
      CALL SECBD2(2,ZML,ETA,NALP)
C
      CALL GIMBL2(2,ZML,ETA,NALP)
C
      CALL GMBLD2(2,ZML,ETA)
C
      CALL SAGIM2(2,ZML,ETA,NALP)
C
      IF(IVISCS.EQ.1) CALL VISCUS(2,ZML,OMEG,OMEGDT,SYDOT,YDOT,MPRIME)
C
      CALL VISCS2(2,ZML,OMEG,OMEGDT,DEPEND,DERIV,VPRIME)
C
      IF(IVISCS.EQ.0) GO TO 170
C
      N=NUMEQS - 4
      DO 160 I=1,3
      LL=N+I
  160  DERIV(LL)=YDOT(I)
       DERIV(NUMEQS)=SYDOT
  170 IF(IDAMP.EQ.0) GO TO 180
       DERIV(10)=DEPEND(11)
       DERIV(11)=ETA(7)
  180 CONTINUE
C
C     LOAD QDERIV ARRAY FOR FLEXIBLE ELEMENTS
C
      IA=IAB-1
      IB=IBB-1
      IC=ICP-1
C
      DO 240 K=1,IK
C
      NKB=NK(K)
      NTW=NKT(K)
C
      NSZ=2*NKB+NTW
C
      IF(NSZ.EQ.0) GO TO 240
C
      DO 220 I=1,NSZ
C
      SUM=0.D0
C
      DO 190 J=1,NSZ
  190 SUM=SUM + DELINV(K,I,J)*RT2(K,J)
C
      SUM2=0.D0
C
      DO 210 L=1,NALP
      SUM1=0.D0
C
      DO 200 J=1,NSZ
  200 SUM1=SUM1 + DELINV(K,I,J)*GAM(K,J,L)
C
  210 SUM2=SUM2 + SUM1*ETA(L)
C
  220 RT(I)=SUM - SUM2
C
C
      NKB2=2*NKB
C
      IF(NKB.EQ.0) GO TO 225
      DO 222 I=1,NKB
      I1=IA+I
      I2=I1+NKB
      I3=IB+I
      I4=I3+NKB
C
       DERIV(I1)=DEPEND(I2)
       DERIV(I2)=RT(I)
       DERIV(I3)=DEPEND(I4)
       DERIV(I4)=RT(I+NKB)
C
  222 CONTINUE
C
      IA=IA+NKB2
      IB=IB+NKB2
C
  225 CONTINUE
C
      IF(NTW.EQ.0) GO TO 240
C
      DO 230 I=1,NTW
C
      I1=IC+I
      I2=I1+NTW
C
       DERIV(I1)=DEPEND(I2)
       DERIV(I2)=RT(I+NKB2)
C
  230 CONTINUE
C
      IC=IC+2*NTW
C
  240 CONTINUE
C
      CALL ACTFLT(2,TIM1,ETA)
C
C
      IF(IOUT.EQ.1) GO TO 260
      WRITE(6,7009) ZLK
      WRITE(6,7010) ZLKP
 7009 FORMAT('0',1X,'ZLK(K)  5 ',1P10E12.4)
 7010 FORMAT('0',1X,'ZLKP(K) 5 ',1P10E12.4)
C
      WRITE(6,8000)
      WRITE(6,9002) DERIV
      WRITE(6,9002) P1,XMB,MPRIME
  260 CONTINUE
C
      RETURN
      END
